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ABSTRACT 


A computer program was developed to find the 
minimum weight design of a stiffened cylindrical shell 
under uniform internal pressure using finite element method. 

The constraints considered were on induced stresses 
displacements and natural frequency. The constant stress 
triangular element was chosen in preference to the more refined 
36 degree of freedom element — namely the doubly curved 
triangle for idealisation of the shell skin, because of the 
limitation imposed on computer storage and time. Axial rod 
elements were used for idealisation of the stiffeners. 

Shell skin thickness and areas of stiffeners are taken as 
design variables. On the basis of parametric study the 
number of design variables was reduced. Optimisation 
results were sought for two different starting points. All 
computations were done on IBM 7044 computer. 



CHAPTER 1 


INTRODUCTION 


With the increase in the size of modern aerospace 
vehicles . , the demand for light weight structures 

also increases. This has made the structural engineers more 
and more conscious of minimum weight design. A structural 
configuration that is widely used in aerospace vehicles is 
the stiffened thin cylindrical shell. 

Stiffened cylinders under bending are representative 
of aircraft fuselages whereas launch vehicles and some space- 
craft can be treated as stiffened cylinders under axial 
compression. 

Leggett and Hopkins - C3j showed that sandwich shell 
construction becomes increasingly more efficient than 
stiffened panel construction as the curvature increases. 
Wittrick |22j investigated the optimum design of a sandwich 
cylinders under compression. He considered various face 
materials and core densities and concluded that there is a 
definite optimum core rigidity for each face material and 
design index (N/d), where N is the load d is the mean diameter 
of the cylinder. Joycee and Mitchell jjlj determined the 
spacing between the stiffeners for minimum weight stiffened 
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cylinders. Gerard and Papirno fio] have considered the 
minimum weight design of ring stiffened and longitudinal 
stiffened cylinders under compression based on orthotropic 
cylinder theory. Richer |12 j treated the problem of optimum 
design of stiffened cylinderical shells subject to multiple 
load conditions. 

The minimum weight design of a cylinder represents 
a constrained non-linear optimisation problem and is usually 
treated by one of the two basic mathematical programming 
methods, namely, the constrained technique or the direct 
search method and the unconstrained or the indirect method. 
The most widely used indirect method is the penalty function 
method. This method consists of solving the constrained 
optimisation problem as a sequence of unconstraint minimi- 
sation problems. 

There are a number of ways in which the unconstr- 
ained minimisation can be accomplished. The Random search 
method, Univarriate method, Powell conjugate direction jl7j , 
Pie tc her and Reeves conjugate gradient method (_6j , Quasi- 

- .... i ' '*} 

Newton method j 3 , or Davidon-P let chef Itowell. variable metric 
method are some of the methods that could be us6d. The 
selection of a particular method depends upon the size (i.e. 
number of variables) of the problem and the difficulty in 
computation of the gradient of the objective function. Por 
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small size problem with a continuous function, Fletcher 
Reeve’s conjugate gradient method or Davidon Fletcher Powell’s 

variable metric method prove to have a faster convergence 

/ 

rate. The variable metric method is numerically more stable. 

The numerical stability of the variable metric 

method depends upon the accuracy with which the minimising 

step length along any particular direction of descent is 

determined. Although the Golden section search method |4j 

has a slow convergence it is the most powerful technique for 

Iter mntixmonm £uaeti»ffi8 

discontinuous functions. r the cubic interpolation 

method is used for one dimensional minimisation for accuracy 
and faster convergence within the variable metric method. 

The evaluation of the penalty function involves the 
evaluation of the constraints. The constraints considered 
in the problem are on the induced stress, natural frequencies 
and the maximum displacement. Thus to evaluate the const- 
raints, a model is required to represent the structure. This 
model or the analytical idealisation must be accurate enough 
and must be suitable for embedment into the iterative design 
analysis loop of the optimisation procedure. 

In this work, the finite element method is used. 

This numerical method is now widely used for structural 
analysis. Within this method are two approaches, namely, 
t the force or the flexibility method and the displacement 
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or the stiffness method of the two, the displacement method 
is easier to use in automated design procedures. Zienkiewich 
j 24 1 and Prezemieniecki [18} describes these methods in great 
detail. 

Two types of elements are used to model the stiffened 
shell. A constant stress element is used for idealisation 
of the skin and social force members for idealisation of the 
stiffeners. The constant stress element was selected after 
a comparative study was made with the doubly curved shallow 
shell triangle element. 

The determination of all the natural frequencies 
is a very time consuming process. Since only the first few 
natural frequencies are needed, it will be a waste of computer 
time to evaluate all the frequencies. The static condensa- 
tion is one method by which the size of the eigen-system can 
be reduced. As the bending degrees of freedom are predominant, 
the inplane degrees of freedom are removed from the eigen- 
system by using the equilibrium equations. Another method by 
which the size of the eigen problem can be reduced is the 
subspace iterative scheme [2]. In this method, a set of 
eigenvectors corresponding to the desired eigenvalues scan 
a subspace. By an Iterative scheme, these set of eigen- 
vectors converge to the required eigenvectors. The attraction 
of this scheme lies in, that in any optimisation stage the 



5 


starting set of eigenvectors are available from the previous 
optimisation cycle. 

In this work a stiffened cylindrical shell clamped 
at both edges and subjected to uniform internal pressure is 
considered. This simplified model is considered due to the 
limitations on computer time and storage on IBM 7044 computer. 
A parametric study is conducted to find the variation of the 
various quantities with respect to the design parameters. 

The results of the parametric study are used in reduc ing 
the number of design variables in optimisation. 

Optimisation results are obtained by considering 
the skin thickness and cross sectional areas of the stiff- 
eners as design variables. 
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CHAPTER 2 . 

FORMULATION AND SOLUTION OP THE PROBLEM 

2.1 EORMULATI OH : 

The objective is to find a minimum weight design of 
a cylinderical shell under uniform pressure. The behavioural 
constraints that must be satisfied are: (1) the stress 
induced must not exceed the yield strength (2) the first few 
natural frequencies must satisfy the specified limitations 
and (3) the maximum displacement must not exceed a specified 
value . 

Mathematically the problem can be stated as follows: 

Minimise P(X) (1) 

0 

> 0 

> 0 ( 2 ) 

> 0 


subject to 

o- y -<r 3 (S) 
-Wk® 


where 
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^(.X) = weight of the structure (objective function) 

sry = yield strength of the material 
(X ) = stress induced in the element 

J 

U? k (X) = k natur a 1 frequency 

Uj^CX) = upper bound on k^* 1 natural frequency 

£t? k (X) = lower bound on k^* 1 natural frequency 

gJZ) = value of the m"^ displacement degree of freedom 
TJl "ttl 

^m = u PP er bound - on the m displacement degree 
of freedom 

X = vector of design variables 

The design variables are taken as the skin thickness 
and the areas of cross-section of stiffeners. 


2.2 INTERIOR PENALTY METHOD OE SOLUTION: 



•X = step length 
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It is in the determination of that the various 
techniques of optimisation differ from one another. One such 
method for unconstrained minimisation is to set 

S q = -VP(I) (4) 

i.e. X q+1 = X q VF(X) (5) 

This can also be written as 

^ = - v»® (6) 

1(0) = x 0 

Thus a minimum could be obtained by solving a set 
of differential equations. The solution X(X”) is called 
"orthogonal trajectory" from the point X Q and it approaches 
a point t/F = 0 as oK, — * oo . 

For unconstrained minimisation equation (6) can be 
modified to incorporate the influence of constraints. One 
approach is to add something to equation (6) which would bend 
the trajectory into the feasible region. This is the under- 
lying principle of the penalty function methods j’8j. 

A constraint divides the n-aimensional design space 
into two parts; the feasible region and the infeasible region. 
Exterior penalty method approaches the minimum from the 
infeasible region, while the interior penalty function 
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method approaches the minimum from the feasible region. In 
this work, the interior penalty function method is used. 

The appealing feature of the interior penalty- 
function method is that given an initial acceptable design, 
it produces an improving sequence of acceptable designs. 
Moreover we approach the constraints in this sequence in such 
a way that they become critical only near the' end. This is 
a desirable feature in an engineering design process because 
instead of taking the optimum design we can choose a sub- 
optimal design but less critical design if we like. Such 
designs have a reserve capacity to absorb overload or abuse. 

Mathematically the interior penalty function can 
be written as 

0(X, r) = F(X) - r £ — (7) 

d=l g-j(x) 

where 

r = penalty parameter 
"till 

g. = 3 constraint function 
D 

m = number of constraints 

and F(X) is to be minimised overall X satisfying g.(X) C 0, 

J 

j = 1,2,..., a. Mote that if r is positive the effect is to. 
add a positive penalty to F(X). This is because at an 
interior point all terms in the sum are negative. As a 
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boundary of the hypersurface is approached, some g. will 

3 

approach zero and the penalty will explode. The penalty 
parameter r is made successively smaller in order to obtain 
the constrain^ minimum of F(X). 

The constraints can be normalised and the minimis- 
ation problem can be restated as follows. 

Minimise F (X ) 

subject to 


ti 


- i 


<T“ 


7 




1 * 3 . 


,U 


0 


< 0 


1 - 


ctsf 

k 


0 


( 8 ) 


e u 


< 0 


- X. 

1 


0 


The algorithm for interior penalty function method 
is given by the following steps. 

(1) Given a starting point X , satisfying all 
g^ (X Q ) < 0 and an initial value of r, minimise 0 to 


obtain X- 


W 
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(2) Check for the convergence of X^„ 

(3) If the convergence criteria are not satisfied 

9 

reduce r as ( r ) new = ^ r ^ 0 ld c w -^ ere c - 

(4) lake X^ as the new starting point for minimi- 
sation and repeat from step 1 onwards. 


2.2.1 DAVIDON -FLETCHER POWELL METHOD : 


The original inventor Davidon j 5j considered this 
method as a terms of a variable method, some consider it as 
a Quasi Newton method and Fletcher and Powell [J \ think it 
to be a conjugate direction method. This method works 
according to the following iterative scheme. 

(1) Start with an initial X q and an initial positive 
definite symmetric matrix H q (for example the identity 
matrix) and set 

S 0 = - [HjVE 0 (9) 

_ _ * _ * 

(2) Compute X . = X + ■-< S where minimises 

q+l q. q. 0. ' 1 

F(X^ + AS ). Check X^ + .j for optimality. 

(3) If X 1+1 is not optimum, compute 


! H 


= ! H_ j + } M ] + 


q+1 ! “ f. qj 


where defining 


LLl = -wd,) 


( 10 ) 
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* 


T 

S S 

_jg g. 

S T Y 

q 2 


(H Y ) (H Y ) 
g a q a 
T 

' Y H Y 

q q q 


T 


(4) Compute S Q+1 = - \ R^ + ~\ V^(X q+1 ) 


from step 2. 


( 11 ) 


and repeat 


It can be shown that if the Hessian matrix [ H q ) is 

positive definite then \ H q+ ^ ] is also positive definite 
* ' _gY _ 

provided X is such that S VF(X .J = 0. The positive 

q q q + 1 

definiteness is preserved in theory only if X is the true 

minimum point, but round off error in other computational 

approximations can cause trouble. Hence in this method care 

must be excercised to ensure that the IhI matrix is not 

updated with data arising from poor approximations used for 

computing \/,' q until S' q V^(^ +1 ) is sufficiently small or 

T —•* 

alternatively skip the update cycle when S q y F(X q+ ,j) ^- s 
large i.e. set [ H q+1 j = j H q ~j and S q+1 = [ H q+1 l V? (£ q+1 ) and 
continue as before. 

Bard [it reports that if the function is badly 
distorted or eccentric or if the function or the variables 
not properly scaled then it is possible that the elements 
of the matrix ( H ~\ and [i 1 would differ from one another 
in order of magnitudes. In finite digit arithmetic this 
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C<2^ 

would imply that one of the matri^ would have no contribution 
to the [H | matrix. This could lead to a non-positive 
definite matrix [h], In such cases the most efficacious 
remedy would be to set |H{ back to [ Hq” 1 and proceed as if 
starting over. This resetting should not be done too often 
otherwise the method cannot be expected to work well. It has 
been recommended that it is a good practice to reset the |H| 
matrix after n-cycles, where n is the number of variables. 

2.2.2 OWE DIMENSIONAL MINIMISATION : 

In the last section it was mentioned that the 

D.E.P. method’s stability is dependent on how accurately the 

* 

minimising step length <X is determined. Since D.E.P. 

requires evaluation of the gradients, it inherently assumes 

that the function is continuous. Hence the slope (first 

derivative with respect to :A) information can be used for 

* 

determination of X • function is approximated by a 

cubic polynomial. 

H (<>.;') = a + bX+ c<* 2 + (13) 

The function and slope values at point A and B 
(see fig. 1) are used to determine the values of a, b, c 


and d. 
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t 


The minimum will he at a point where 
which yields 


dHQQ _ n 
d>( “ u ’ 


h + 2c <+ 3d* 2 = 0 


(H) 


- c + (c 2 - 3hd) s 

or X" = =-33 (15) 

2 

c\ H 

As — p has to he positive at the minimum, the positive sign 

dor 

has to be taken in Eq. (15). The constants a, h, c and d can 

he evaluated in terms of F(A), P(B) and P 1 (A) and F'(B), 

where F(A) , F 1 (A) and F(B), P’ (B) are the function and slope 

values at points A and B respectively. To ensure that 
2 

c - 3bd is greater than zero, we have to make sure that 
F’(A) 0 and P 1 (B) 0. Substituting the values of a, h, 

c and d in equation (15) and rearranging the terms, ■.< can 
be expressed as: 


* 


A . P» (A) +Z_+Q 
A + P» (A) + P' (B) + 2Z 


(B - A) 


where 


Z* = 3 jA . )f-^ i, g . )J, + pt (a) + p> (B) (16) 

Q = (Z 2 - P' (A) p« (B))* . 
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CHAPTER 3 

ANALYSIS OP THE STRUCTURE 


To evaluate the stress, displacement and frequency 
constraints, the analysis of the structure is required. To 
analyse the structure, the finite element method is used. 
This method consists of representing the structure hy an 
assembly of elements as shown in Eig. 2. 

Two types of elements are used to represent the 
shell in this work. The triangular elements are used to 
represent the skin and the axial force members are used to 
represent the stiffeners. 

To represent any general shell structure accurately 
the doubly curved triangular elements capable of taking 
bending and inplane loads have to be used. However for 
simplicity the flat triangular elements, capable of taking 
constant stresses, are used in this work. 

A comparative study is made here to find . the 


efficiencies of a doubly curved shallow shell element [15,16] 
and the constant stress flat triangular element with respect 
to accuracy, rate of convergence and storage required. 

Before presenting the numerical results of the comparative 
study, the derivation of the stiffness and the mass matrices 
of the two element is outlined. 



16 


The elemental stiffness and mass matrices are 
derived from the local assumed displacement polynomial. Thus 
the finite element method may be thought of as an assumed 
mode analysis in which the assumed mode instead of having a 
uniform definition for each element has a sort of patch work 
of assumed modes. The assumed polynomial has a number of 
unknown parameters which are evaluated in terms of the nodal 
degrees of freedom of the element. By using the strain- 
displacement and stress-strain relations of linear elasticity 
the strain and kinetic energies of the element are expressed 
in terms of the unknown nodal degrees of freedom. 

By expressing the strain and the kinetic energies 
as quadratic forms in the nodal displacement and nodal 
velocity vectors respectively, the element stiffness and mass 
matrices can be identified. The element stiffness and mass 
matrices are usually determined with respect to the local 
co-ordinate system. These matrices are then transformed to 
correspond to a common global co-ordinate system before 
assembling the structural stiffness and mass matrices. The 
assembled stiffness and mass matrices are used to state the 
equilibrium and the eigenvalue problems in the standard 


manner. 
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3.1 DOUBLY CURVED TR1AUGLE (D.C.T.) ; 

The geometry of an arbitrary triangular shallow 
shell element is shown in Pig. 3. 

Here u, v and w are the displacements in the local 
co-ordinate system y , - , Z and X, Y, Z represents the global 
co-ordinate system. The dimensions a, b and c and the 
rotation angle © are defined in terms of the global co- 
ordinates of the vertices. f (2, q) defines the shape of the 
shell about the base plane y_- v ; . 


3.1.1 STIPPNESS MATRIX; 

Prom the shallow shell theory of Uovozhilov I 1 4 j , 
the strain energy is given by 


V 


r 


Bt 


fr r 2 

1 } U- 


2(1 ~ V ) J) 1 


~ 2 [( j\-j+ V j + 0 -I 1 ) ^ y\ 


+ v 2 + 2 yu ,Vy^ + ^(1 - v) (u y . + v- £ )‘ 


(u n + } w + [ (y -£^ + y hn + 2 V ^ W 


2(1 - V) . 


c 2 ■> 2 


' in 


w + 


n.n 

t 2 r _2 


72 I w t. t + w ’“. n + 2 ^ w n w nn 


+ 2(1 -V0w| n ]}d.^dyj 


( 17 ) 


where the integration is to be done over the base plane 
1 f 2 1 3 * . The subscripts of u, v, w and f* represents the 
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derivatives with respect to tfand v<. For example 


u. 


an 


w. 


2 

Tw 


etc . 


^ .yi. a- 

The function S (^,\|) is assumed to be of quadratic 


form: 


T (£.»*}) = C 1 + °2*~ + c 3 > ’\ + c 4 + c 5't,.^ + c 6 rj 08) 

This implies that the shell element has constant 
curvatures and this is consistent with the shallow shell 
theory. 

A quaiyfcic polynomial is assumed for the normal 
displacement w and a cubic polynomial for the tangential 
displacement. Thus the tangential and the normal displace- 
ment yield the same energy convergence rate n~^ where n is 
the number of elements. 

2 2 3 

u( K|) = a^ + a 2 £_+ Y 1 "*■ a^ £_ a g t ^ Q-g t 

2 2 3 

+ a 8l ^ + a 9 ^ ^ + a 'io , '\ 09a) 

2 2 

v(?.,iq.) = + a 12 £ + a^.;j+ + + a i6’ r ’| 

+ a 1? £ > + a 18l _ 2 ^+ a 1g £^ 2 + a 2 Q t ^ 3 


($9b) 
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r (Z>*f) 


a 2l + a 22 i . + a 23 > ] + a 24^ + a 25£ T l + a 26 y f + a 27?_ 5 


+ a 28sL>) + a 29^ Y 'i + a 30^ + a 3ll + a 32‘l Y \ 

+ s 33i*f + a 34 ~ r f + a 35 'f + a 36l + a 37^ r f 


+ a 38l 2? ? + a 39^n 4 + a 40^ 


( 19 c) 


The nodal displacement degree of freedom for 
u ( £. ,V f ) and v(£,)i) are taken as u and v and their first deri- 
vatives with respect to f and at each of the corners plus 
u and v at the centroid of the element. This gives a total 
of 20 degrees of freedom for the element consistent with the 
20 free parameters in equation (19a) and (I9h). The nodal 
displacements for w(7_ r 4^ are taken as w and its first and 
second derivatives with respect to ~p_and \~>at each of the 
three corners. These along with the constraint that the 
normal must be cubic in variation along the edges, yield 
sufficient information for evaluating the 20 free parameters 
of (19c). All generalised displacement are assembled into 
a 38 column vector J in local co-ordinate system £ , 
as 
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The subscripts 1, 2, 3 and c denotes the comers 
1 , 2 and 3 and the centroid of the element respectively. 

The coefficients of equations (19) are assembled into a 

column vector A l as 

} ~ ? ®"2 J • • • • j (21 ) 

Combining equations (19) to (21) one obtains 

A W 1 } 

f 0 [• = .\t) \a\ ( 22 ) 

. ^ 

V o j 

where the 40 x 40 T matrix is given in Appendix A. The 

34 31 2 2 

determinant of T has a value -640^ (a + by (a + c ) 

(b + c ) /729 which is non-zero for all practical problems. 
Hence the equation (22) may be inverted to obtain 

[w 1 ) 

fAl = [t" 1 ] jo l = 1^] 4WJ (23) 

.1 i . 

0 I 

where the 40 x 38 matrix \ T. ] consists of the first 38 columns 

of It- 1 :. ' ' 

The displacement functions of equations (19) can 


be rewritten as 
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10 m. . n. 

u = a. ^ ^ 

0=1 3 


20 p . q . 

Y = a. X 3 v ; 3 (24) 

3=11 3 


40 r. s . 

w = y a . ~A. ^ ' i ^ 

3=21 ^ 


By substituting the displacement function of Eq. 
(24) into the strain energy expression (17), and performing 
the integration, the strain energy can be expressed in 
matrix form as 




(25) 


where the elements of the 40 x 40 matrix [k] are given in 
Appendix A, Combining equations (23) and (25), the strain 
energy can be obtained in terms of { w as 


U = — p- hfljE ])w\ (26) 

2('t - V ) 1 •' L V. t- 1 1 

where i K-. | = [ic ] is the 38 x 38 stiffness matrix 

in terms of the generalised displacements relative to local 
co-ordinate system using the relation. The stiffness matrix 
relative to the global co-ordinate system can be obtained 
using the relation 
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= [e1{W 2 | . (27) 

where 

{ W 2) = (u 1> 'M’ u y1’ T 1’ T x1’ V y1 ’ w 1> m x1’ W y1 ’ 

w xx1 J w xy1 5 w yy1 J U 2 J ***’ U 3 J *'*’ u c* v c^ 

(28) 

where ) Wg [, is the generalised displacement vector relative 
to the global co-ordinates and I R j is the rotation matrix 

i , 

given in Appendix A. Equation (26) and (27) yield 


u 


Et 

2(1 ~Y' 2 ) 



f L K 


2 j 



where the matrix 


(29) 



(30) 


is the 38 x 38 stiffness matrix relative to the global 
co-ordinate system. 


3.1.2 CONSISTENT LOAD VECTOR: 

The consistent load vector is obtained by calcul- 
ating the virtual work done by the applied loads. Q^Cl,^), 
Q v (Yj-'i) arj d QyjCl.jh) in the u, v and w directions respectively. 
The transformed load vector becomes 




(3D 
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where "the entries of |Q"jj s^e given by 


K 

m. 

X 1 

ru 

V) 1 

v 

dn 

d V; 

i 

o 

Vi 

S.K 


•H 

dl. 

d n 

10 

< i < 20 

!K 

r. 
X 1 

Vt S i 

dx 

d y \ 

20 

< i r 40 


(32) 


In the case of a uniform pressure load of intensity 



f 0 
\ 

J 

■y 





i < 20 
i > 20 


( 33 ) 


The centroidal displacement degrees u and v are 
condensed out as they lie inside the element. This reduces 
the number of degree of freedom in each element and also 
simplifies the assembly procedure. The reduction is carried 
out as follows. 

The elemental equilibrium equations can be written 


in partitioned form as 



i 

i 

i 


\ 






( 34 ) 
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where j W q ^ is a 2 component vector (u c , v Q ) an d j ^ 0 \ 
contains the first 36 component of equation (28). Equation 
(34) can be written as two separate equations and W can 
be determined to obtain 

[K']4Wi = (35) 

where 

[k| = f K 1 - [ K 1 1 " K -1 1 j~K ~\ T ( 36 ) 

1 L o J i ocj <■ c Ji L ocj w ' 

is the reduced 36 x 36 stiffness matrix and 

M - [ p °!r ■ W K 1 ] l p c] (37) 

is the reduced consistent load vector. 

3.1.3 MASS MATRIX: 

The kinetic energy of the shell element can be 
expressed as 

T = -S-| j \(u 2 + v 2 + w 2 )ds. dy^ (38) 

As the shell vibrations involve predominantly 
normal motion, within the domain of small deformation theory 
the tangential inertia can be neglected. Assuming harmonic 
time dependence, Eq. (38) can be written as 

T = j I w 2 d * dr^ ( 39 ) 
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where co is the circular frequency. Substitution of Eq. (24) 
into Eq. (39) yields 

= " L ir~~ i A \ T L m l \ A \ (40) 

F 0\ + nr, n ± + n..) i >20, j > 20 

(41) 

0 i < 20, 5 < 20 

Following the same transformation procedure 
described for the stiffness matrix, the elemental mass ■ 
matrix in global co-ordinate system can be obtained as 

l M ] = [ r J T [ T i] T [*] I T l] [ E ] (42) 

3.2 CONSTANT STRESS FLAT TRIANGULAR ELEMENT. (C . S .1 . ) : 

This element has a constant stress field within 
the element. By Hook’s law therefore the strains are also 
constant. From strain displacement relationship we get a 
linearly varying displacement field. It is also assumed 
that the element can only undergo inplane displacement as 
shown in Fig. 4. 

3.2.1 STIFFNESS MATRIX ; 

The local displacement variation can be assumed as 


T 


where 


m. 


13 
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U.(aL > "V|) — a-| + Eg > + a^ Vj 

v(fc, V^) = + a^~5_. + ag r«. 


(43a) 

■(43b) 


The constants in Eq. (43) can be evaluated in terms 
of the nodal displacement degree of freedom and the nodal 
co-ordinates and Eq. (44) can be written in matrix form as 


u "} 


. Y 


v. = ! a] 


i* j r i 1 


(44a) 


where the elements of the matrix j A j is given in Appendix B 
and W ^ is the vector of generalised displacement relative 
to the local co-ordinates 


\ W 1 l = j Vy, "V-j ? u 2 , v 2 , u 3 , v 3 ] 


. '' s 


J 

The strains can be obtained by differentiating 
equation (44) to yield 




l T = [b] 


(45) 


T . 


where the matrix > b j is given in Appendix B and • £ [ is 


i, ) 


the strain vector i 


i -a- .... Li.. . ■ ' " 'V* 

v <: *v. •* ^ * ( V; *. v. 7?; y* ; 

v i «c * l i 

By Hook's lay/ the stresses can be written as 


V. ' ) 


i c ] 


L 


(46) 
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where the 3x3 matrix [c] is given in Appendix B and 
is the stress vector 

The strain energy of the element is given by 

U = | hr-l d *H (47) 

By substituting equations (45) and (46) in equation 
(47 ) , the strain energy can be written as 

u = | ( [) w i“i T r b i T ! c ] w \ w i i d n (48) 

As the matrices \b] and \o \ and the vector) W. t 
are constants, they can be taken out of the integral and the 
integration can be performed to obtain 

V - 5{ w i] t l K iH w n (49) 

where the 6x6 matrix 

' K n l = /\t [bj T [cj [bj (50) 

is the stiffness matrix relative to the local system and A 
is the area of the triangle. 

This matrix is now transformed to the global 
co-ordinate system using 

(V = [*J ?* 2 \ ' (51) 

where the 9x6 matrix f R ] is given in Appendix B and 
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| W 2l - | U 1> Y l> w r , u 2 , u 3 ...| T ( 52 ) 

is a vector of generalised displacements relative to the 
global co-ordinate system. 

Substitution of Eq. (51 ) in Eq. ( 49 ) yields 

f K 2 ] = | R'p | IC 1 < | R | - (53) 

where the 9x9 matrix 1 Kg j is the stiffness matrix relative 

* -» 

to the global system. 


3.2.2 MASS MATRIX : 


The kinetic energy of the element is given by 


® = \ \ \ (u 2 + v 2 ) dR d r ( 


(54) 


Assuming a harmonic time dependence and substitution 
of Eq. (43) into Eq. ( 54 ), one can express the kinetic 
energy as 



(55) 


where the 6 x 6 matrix j m ) is the mass matrix relative to 
the local co-ordinate system and is given in Appendix B, 
and u : is the circular frequency. 

following the same transformation procedure as 
given for the stiffness matrix, the 9x9 mass matrix \ M ] 
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relative to the global co-ordinate system can be obtained 
as 

re] = jRftm] [E] (56) 

3.3 EVALUATION OF C0FSTBACT1S ; 

Once the elemental stiffness and mass matrices are 
available, the stiffness and mass matrices of the whole 
structure can be obtained ^18,24 j. After introducing the 
boundary conditions, the displacement and the natural 
frequencies can be obtained respectively by solving the two 
system of equations: 

[k]\w! = (57) 

(|_K] - 68^ [Mj) \ W] = 0 (58) 

where [ Kl and 1 M 1 are the structural stiffness and mass 
matrices respectively, \ W \ is the displacement vector and 
j P ( the external load vector and ix> is the circular, frequency. 

Once the displacements are known the stresses can 
be calculated. The principal stresses are found and 
Von-misses yield criterion is used to evaluate the stress 
constraint. 

To evaluate the natural frequencies, the subspace 
iteration scheme described in the following section, is used. 
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3.4 SUBSPACE ITERATIVE SCHEME f 2? ; 

The objective is to solve for p— lowest eigenvalues 
and the associated eigenvectors satisfying the relation 

= MS* 1 }! a 2 ] (59) 

where j 0 stores a set of p-eigenvector corresponding to the 
frequencies [ A 2 ]is a diagonal matrix j A 2 "j= ^{j?;- 2 J 
where are the rotational frequencies. 

The essential idea is to iterate simultaneously 
with p-linearly independent vectors which initially span the 
starting subspace S Q , until S' is spanned. The required 
eigenvectors are then computed without further iterations. 

The total number of required iterations depends upon how 
close S Q is to S . But the effectiveness of the algorithm 
lies in that it is easier to establish p-dimensional 
starting subspace S Q close to than to find p-vectors, 
each of which is close to the required eigenvectors. Also 
convergence of subspace is all that is required and not of 
individual iteration vectors to the eigenvectors. 

The algorithm is given by the following steps; 

(1) Let | store the set of p-eigenvectors 
corresponding to the natural frequencies desired. 

(2) Solve the set of simultaneous equation 

AllVj - K AA 


(60) 
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(3) Find the 


projection of the operations f K ~\ and 


[ffi.jonto Sj 


FkL = )x/i T f 


*k “ \ ~ k> 




\ M 1 = i X, l ; M 1 ) X. I 
•- k i kj l J i k ; 

(4) Check for convergence i.e. J ffi ■ ' i = F 1 3 i.e. 


(61) 


( 62 ) 


identity matrix. 


(5) Solve for the eigensystem of the projected 


operators 


l K SM < 


(63) 


(6) Find an improved approximation to the eigen- 


vectors 


X k - 


and return to step 2. 


Then provided S is not orthogonal to one of the 


required eigenvectors 


0 as k gd is obtained. 


Assuming that S^. is close to , the convergence 
“til 2 2 

rate of the i column in X^ is lo -j . Although this 

is an asymptotic convergence rate, it indicates that the 
lowest eigenvalue converges at the fastest rate. As the 
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vectors span S^., they become more and more parallel and 
hence form poorer and poorer basis. One way to overcome 
this difficulty is to generate orthogonal bases in the sub- 
space S k . Another way to accelerate the convergence is to 
choose q eigenvectors with q > p and to stop the iteration 
process whenever p of the eigenvectors have converged. The 
originator of this method Bathe and Wilson '^2 recommended 
a value of q = min (2p, p+8). 

To find the starting subspace, the following is the 
recommendation of Bathe and Wilson [ 2 \. The first column of 
MX q (Eq. (60)) is set equal to the diagonal elements of j_Mj. 
This assures that all degree of freedom corresponding to 
the natural frequencies are excited in order not to miss a 
mode. The other columns are not vectors with +1 at the 
co-ordinates. with the largest ratio. 

3.5 C OMPARATIVE STUDY ; 

A comparative study between constant stress trian- 
gular element (C.S.T.) and doubly curved triangular element 
was made using a simply supported spherical cap shown in 
Eig. 5, as a test problem. The cap was subjected to uniform 
pressure. As the structure and the loads on it are symmetric 
in nature only one quarter of the cap is taken. The 
idealisation of the structure using D.C.T. and C.S.T. are 
shown in Eig. 6(a) and 6(b) respectively. 
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Sixty three degrees of freedom were taken for the 
idealisation with D.C.T. and seventyfive degrees of freedom 
for the idealisation with C.S.T. The non-dimensional 
parameters (Ety/p Q El) and (Etw/p^R^) are plotted along the 
edge AB and BC (see Eig. 5) respectively in .jtig- 7 and 8 
From the plot ; . it is apparent that even with lesser degree 
of freedom. D.C.T. gives more accurate results than C.S.T. 

The execution time on IBM 7044 computer using 
D.C.T. for idealisation for one static analysis and one 
dynamic analysis for the two lowest frequencies was found 
to he 5 minutes 3 seconds, while for C.S.T. it was found to 
be 51 seconds. Mow on an average, an optimisation problem 
requires nearly 600 function evaluation (Ref. ) 12l). Thus 
only time consideration prohibits the use of D.C.T. for 
idealisation of the skin of the shell. 

The other problem with D.C.T. is that of storage. 
D.C.T. has twelve degrees of freedom per node, while C.S.T. 
has three degrees of freedom per node. Thus for every node 
of D.C.T., four nodes of C.S.T. can be taken. Mow in an 
idealisation of the stiffened shell the most convenient 
idealisation is that of placing a node at the intersection 
of a circumferential and longitudinal stiffener. Ref. \l2j 
reports that near the optimum shell design, the longitu- 
dinal stiffeners are separated by two or three degrees from 
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each other, Thus the number of nodes can be very high and 
hence the storage problem in case of D.C.T. will be very 
severe. 

Thus even though D.C.T. yields more accurate results 
than C.S.T., the time and storage consideration shows that 
C.S.T. should be used for idealisation of the shell skin. 

3.6 CONVERGENCE STUDY OF C.S.T. ■ 

Dor the spherical cap problem a number of ideal- 
isation were made as shown in Pig. 9. The variation of the 
non-dimensional displacements (EtV/Rl) and (EtW/R ) (along 
the edge AB and BC respectively) with the number of elements 
are plotted in Pig. 10 and 11 respectively. Por the ideal- 
isation with elements 4, 6 and 8, the results obtained are 
found to have no relation with the solution and hence are 
not plotted. The plot in Pig. 10 and 11 show that with the 
increase in number of elements, the results tendjb towards 
the exact solution. In Pig.. 12, the variation of the two 
lowest natural frequencies with the number of elements are 
plotted. Prom Pig. 10, 11 and 12 it can be seen that the 
structure becomes less stiffer with the increase in number 
of elements. 

The number of elements is inversely proportional 
to the step size (length for one dimensional elements, area 
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for "two dimensional element; and volume for three dimensional 
element). As the number of elements increases the step size 
decreases. The truncation error (in this case the difference 
between the actual mode shape and the mode shape obtained 
from the patch work of the approximated mode shape of 
element) decreases with the decrease in step size. But the 
round off error caused by the finite digit arithmetic of 
the computer, increases with the decrease in step size. To 
obtain the optimum number of elements which would yield 
the minimum error, the results must be plotted against the 
number of elements and that point should be tuken where the 
results do not change appreciably with the change in number 
of elements. Unfortunately this could not be done as the 
storage limitation of IBM 7044 computer is 32 K bits* Hence 
the absolute values of displacements, stresses and 
frequencies must be taken with caution. Yet behaviour of 
these quantities at any node with respect to other nodes 
can be taken. Also the behaviour of these quantities with 
respect to design parameters can be taken. Bor in both the 
cases the effect of the error would be uniform. 
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CHAPTER 4 

RESULTS AND DISCUSSIONS 


In any optimisation problem, the time consuming 
operation is the function evaluation. A function evalu- 
ation requires the analysis of the structure. To save 
computer time, the displacements and frequencies are 
linearised over a small domain of the design variables, 
thereby avoiding the regular function evaluation until the 
design variables change by a certain percentage. 


4.1 DISPLACEMENTS AND NATURAL FREQUENCY GRADIENTS: 


Assuming that the load vector is independent of 
the design variables and the mode shapes change very little 
with a small change in the design variables, Eqs. (57) and 
(58) can be differentiated with respect to the design 
variables to yield. 


dX v J cjx 


= 0 


( 


;aw 


2 tx4l'- CO ? TT—) w \ = 
<5 x ^ v t. i 


<?x ’ Sx 

r -i i 

Premultiplying Eqs. (64) by [ _K J and Eqs. 
and rearranging the terms one obtains. 


(64) 

0 (65) 

(65) by{w} T 
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dX 



(66) 


> f k] 


->x 


W'i T TTT- ] W • -U) 2 | w'f jw l 

dl ' 1 ' H X ' ' 

2u/:W> T H1 !w( 


(67) 


Thus the evaluation of the displacement and 

f K] 

frequency gradients require the evaluation of and 

3L m j m , , J sl k j ;U®i 3 s 

—3—. The evaluation of — — and — — requires the 

d X c X ,-) X 

assembly of the stiffness'-;, and mass matrices at the per- 
turbed point if a finite difference scheme is used. Thus 
we do not have to solve an eigenvalue system or a set of 
simultaneous equations to evaluate the natural frequency 
and displacements- gradients. Once these gradients are 
known, they are assumed constant over a small domain of 
the design variables and the displacement and frequencies 
are calculated directly using the relations. 


\ W l 


mew 


n > 

= Hold + TV 


o. 


( 68 ) 
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n rHO 
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(69) 
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4.2 CYLINDRICAL SHELL PROBLEM: 

The problem considered is the minimum weight design 
of a thin cylindrical shell clamped at both ends and subje- 
cted to a uniform pressure as shown in Pig. 2. The skin of 
the shell is idealized by constant-stress triangular elements 
and the stiffeners byaxial rod elements as shown in Pig. 1. 
The design variable s are skin thickness and the cross 
sectional areas of the stiffeners. The constraints are 
placed on the induced stresses, natural frequencies and 
displacements. 

As the structure and loads are symmetric in nature 
only one quarter of the shell' is considered*. The objective 
is to minimise the function 

M M+N 

P(X) = hlPRlX(l) + fL 5. X(i+1) + fjTR ■£. X(j + 1)/ 

i=1 j=M+1 

fTlRL (70) 

where 

X( 1 ) = skin thickness 

X(i+1) for i = 1,17 = longitudinal stiffener areas 
X(j+l) for j = M+1 , M+F are circumferential stiffner 
areas 

M = number of longitudinal stiffeners 

H = nu m ber of circumferential stiffeners 

* At X = 0 and L, u, v and w are zero? At Z = 0, ¥{is zero? 
and at Y =0, w is zero. 
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The upper bound on the natural frequencies and 
displacements are taken as thirty percent more than the 
values calculated at the starting point. 

4. 2.1 PARAMETRIC STUDY? 

A parametric study is conducted for a stiffened 
cylindrical shell with five longitudinal and four circum- 
ferential stiffeners in a quarter shell. In Figures 13 to 
21 are plotted the displacements and mode shapes at the 
reference point. The reference point is taken as 
X(l) = 0.004 inches 

X(2) = 0.010 inches 

X(3) = 0.010 .inches 

X(4) = 0.010 inches 

X(5) = 0.040 inches 

X(6) = 0.040 inches 

E = 10 7 psi 

) = 0.101 lbs/ cubic inch 

'/) = 0.333 

R =150 inches 

1 =60 inches 

Figure 13 is a plot of u, v and w vs. 9 at 


x = 12" for a uniform pressure load of 3". psjt. ; ;The dis- 
placements are sinusodial in nature. The variation of 
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displacements u has a wavelength of radians while the 
wavelength for displacements v and w is 2 TT radians . The 
displacement v and w vary with a phase difference of TT/2. 

This phase difference is in accordance 'with the boundary 
conditions imposed due to symmetric nature of the problem. 

The boundary conditioned imposed are that the displacements 
v and w become alternatively zero after each angular rotation 
of /// 2. Hence the phase difference between v and w must 
be an integer multiple of TT/2 . 

Figure 14 is a plot of u, v and w vs. x at 9 = 0° 
and 22.5°. These plots show that the nature of variation 
of u, v and w vs. 9 is sinusodial in nature for any given 
value of x, but the magnitude of the amplitude varies. The 
variation of the amplitude is symmetric about the mid section, 
that is x = 30 inches. This is also expected in this case 
as the structure and the loading are symmetric about the 
mid sections. 

Figures 15 to 18 are the plots for the first mode 
of vibration. The natural frequency for the first mode of 
vibration was found to be 89.7 cycles per second. The 
displacements v and w arc out of phase by 77/2 radians. 

Figures 19 to 21 are the plots for the second 
mode of vibration. The natural frequency for the second 
mode of vibration was found to be 92.3 cycles per second. 
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The results show that the phase difference 
between the displacements v and w is 377/2 radians. 

Figures 22 to 25 show the plots of the displacements 
and frequency versus the design variables . The numerical 
values are given in Table 3. 

From Figure 22 it is seen that with the variation 
of skin thickness, X(1), by 36 percent the variation in 
the displacements u, v and w are 32, 18.2 and 18.2 percent 
respectively. 

Figure 23 shows that with a variation of area of 
longitudinal stiffeners, X(2), by 45 percent the variation 
of the displacement u is of the order of 2 percent. The 
displacements v and w are practically unaffected. The 
variation of the other longitudinal stiffener areas do not 
affect the results appreciably. This shows that for a shell 
under uniform pressure, the longitudinal stiffeners have 
a negligible effect. 

Figure 24 indicates that with a variation of the 
area of circumferential stiffeners, X(5), by 36 percent, 
the displacements u, v and w vary by 1.7, 2.7 and 2.7 
percent respectively. Figure 25 shows that with variation 
of X(6) by 36 percent the variation in the displacements 
u, v and w aro 13.0, 11.8 and 11.8 percent respectively. 
These results show that the circumferential stiffeners 
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located farther away from the ends of a cylinder affect 
the results to a greater extent compared to the one 
located near the ends. The variation in the two lowest 
natural frequencies is negligible. 

4.2.2 Minimum weight design; 

A number of numerical difficulties were encoun- 
tered in the solution of the optimisation problem. Initially 
the problem was worked with six design variables. The 

first difficulty that arose was due to the approximate 

\ . 

evaluation of the displacements and natural frequencies 
using equations (68) and (69). Towards the end of one 
dimensional minimisation the program got bogged down as 
the computer was evaluating the function at two points only. 
These points were separated by a distance equal to the 
domain over which the displacements and frequencies were 
approximated. This difficulty was overcome by reducing the 
domain of approximation after every five iterations. 

Inspito of this the convergence was so slow that not even 
one unconstrained minimisation could be done in 30 minutes 
of computer time. To accelerate the process, the number 
of design variables was reduced to three. This was done 
by assuming the area of cross section of all longitudinal 
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stiffeners to be same. A similar assumption v/as made for 
the circumferential stiffeners. This increased the conver- 
gence rate slightly but a new difficulty arose. 

Though the penalty function was decreasing the 
objective function was increasing. Two steps were taken to 
overcome this. The initial value of the penalty parameter 
was decreased and the variables were non-dimensionalized 
with respect to the length of the starting design vector. 
After an initial rise in objective function, the value 
started decreasing and at the end of 30 minutes a 20 per- 
cent reduction in the objective was found and all displa- 
cement constraints were activated. The starting and the 
final values of optimisation are given in Table 4.1. In 
Table 4.2 the initial values and the final results of 
optimisation, obtained with a different starting point are 
given. The upper bounds on the displacements and natural 
frequencies are taken as 30 percent more than the values 
corresponding to the starting point. In both these examples, 
fourteen constraints were considered. The first seven 
constraints were on induced stresses. The next four were 
on the natural frequencies. g(12), g(13) an& g(l4) are 
constraints on the displacement u, v and w respectively. 



CHAPTER 5 


CONCLUSIONS AID EEC QMMEND AT ION S 

5.1 CONCLUSIONS: 

The results presented in this work show that 

whenever finite element method is used the optimum number 

of elements which would yield a minimum error, should be 

determined. In case this is not done the error in the 

analysis results may be as high .50 to 60 percent. 

wriati** if mmw m Mm i l mm ** mMiimm 

aflat* m •***»**>" ****** *** **^ , * i * # 

s»» » m*m mmm ?*•»•***. 

The circumferential stiffeners located further 
away from the ends of the cylinder have a greater effect 
on displacements and stresses than the ones located near 
the ends. 

5 . 2 RECOMMENDATIONS s 

In this work the buckling and flutter 
constraints were not included. To include these const- 
raints the spacing of stiffeners must also be taken as 
variables. If an optimum is sought., donsidering these 
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constraints , then near the optimum, the spacing between 
the longitudinal stiffeners will be of the order of two 
to three degrees (Ref. M2]), Since the most convenient 
idealisation involves the placing of nodes at the inter- 
section of longitudinal and circumferential stiffeners, 
this will require a very large number of elements. This 
would require a large amount of computer storage and time 
for tho analysis of the structure. This is the reason why a 
a simplified model was chosen in the present analysis. 

This difficulty can be surmounted by taking the following 
steps. Instead of using two types of isotropic elements 
for idealisation, use an orthotropic shell element. This 
element could be obtained by averaging the elastic 
properties of stiffeners over the stiffener spacings. 

This would help in reducing the number of elements for 
idealisation. By using sparse matrix method instead of 
band width method, the severity of storage and computer 
time problem could be reduced. 

To save the computer time further the optimisation 
technique presented in Ref. j 23] is recommended. In this 
scheme tho optimum is sought in two stages, first, by 
a proper grouping of the design variables, the number of 
parameters that optimise the weight is minimised. On 
the basis of this, a mathematical search technique is 
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employed, and design charts are prepared which clearly 

show that the effect of these few parameters on the weight 
of the shell. 












Fig, 4 Constant Stress Triangle 
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Fig. 8 Variation Of (Etw /p 0 R z ) Along Edge BC 
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Number Of Elements ► 

¥iq.\Z Variation Of Two Lowest Natural Frequency With 
Respect To Number Of Elements. 
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TABLE 1 . COMPARATIVE STUDY 


i 

JD.O.E.* 

CO 1 

CO 2 

Computer time 

! 

! 

cycles/sec . 

L 

cycles/ sec . 


C.S.T. 75 

254 

280 

51 sec. 

D.C.T. 63 

232 

239 

303 sec. 


i 


! 

Changes | 

(EtV/p Q RL) 

» t 

{Changes ] 

(EtW/p R ; 

y = -2 : 

1 

C.S.T. | 

y 

D.C.T. 

tin y {- 
]x = 0 ! 

? i 

c.s.t. ; 

t 

D.C.T. 

- 2.0 

0 

0 

-2.0 

0 

0 

-1 .6 

-0.1613 

- 

-1 .6 

-0.9044 

- 

-1.2 

-0.2618 

- 

-1.2 

-1 .0091 

- 

-1 .0 

- 

-0.3258 

-1.0 

- 

-1.1609 

- 0.8 

-0.3225 

- 

-0.8 

-1.0286 

- 

- 0.4 

-0.3619 

- 

-0.4 

-1 .0200 

- 

0.0 

-0.3850 

-0.4080 

0.0 

-0.8004 

i 

. 

O 


* Degree of freedom 
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TABLE 2. CONVERGENCE STUDY 


Table 2.1 (EtV/p Q RL) variation along edge AB 




50 


- 2.000 
- 1.600 
- 1 .500 
- 1.333 
- 1 .200 
- 1 .000 
- 0.800 
- 0.666 
- 0.500 
- 0.4000 
- 0.009 


0 0 0 0 0 0 

_ -0.1613 

-0.1469 - -0.1773 -0.1770 

-0.1939 - 

-0.2618 

-0.2392 - -0.2796 

-0.3225 

-0.2913 - 

-0.2621 - -0.3423 -0.3360 

_ - - - -0.3619 

-0.3249 -0.3109 -0.3362 -0.3584 -0.3671 -0.3850 
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Table 2.2 (EtW/p Q R 2 ) variation along edge 


Jo . elements 
y-coordina^ 

- 2.000 
- 1 .600 
-1.500 
-1.333 
- 1.200 
- 1.000 
-0.800 
- 0.666 
- 0.500 
- 0.400 
- 0.000 


0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 

- - 0.9044 

- 0.7276 - - 1,1032 0.9138 

- 0.9307 - 

- 1.0091 

- 0.8255 - - 1.0190 

- - 1.0286 

- 1.0142 - 

- 1.0989 - - 1.1458 - 1.0225 

----- - 1.0200 

- 1.0793 - 0.8754 - 0.8113 - 1.1643 - 0.8081 - 0.8004 


Table 2.3 O ^ and&? 2 variation 


No. elements J 10 


10 

i 

{ 14 

t 

! 

j is 

f 

? 

j 22 

i 

* 

32 

■ . 

50 

266 

266 

264 

250 

257 

254 

345 

345 

331 

305 

297 

280 
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TABLE 3. PARAMETRIC STUDY* 


1 

Design {Percentage 
vari- {change in 
able {design 

{ variable 

1 

“1 

1 

JuxlO 

i 

i 

i . » 

{ inches 

i 

! 

1 ! 

vxlO \ 

1 

1 

! 

inches | 

1 

wzIO - "^ { 

t 

t 

r 

inches 1 , 

0 1 

cycles / { 
sec . ! 

^2 

cycles/ 

sec. 

X(1) 

0 

-1.81 

-6.64 

-6.64 

89.7 

92.3 


-20 

-2.06 

-7.70 

-7.70 

89.6 

91.6 


-15 

-1.99 

-7.40 

-7.40 

89,7 

91.8 


-10 

-1.92 

-7.13 

-7.13 

89.7 

92.0 


- 5 

-1.86 

-6.88 

-6.88 

89.7 

92.2 


+ 5 

-1.75 

-6.42 

-6.42 

89.8 

92.4 


+ 10 

-1.70 

-6.21 

-6.21 

89.8 

92.5 


+15 

-1.65 

-6.02 

-6.02 

89.8 

92.7 


+20 

-1.61 

-5.03 

-5.03 

89.8 

92.8 


+25 

-1.56 

-5.66 

-5.66 

89.9 

92.8 

X(5) 

0 

-1.81 

-6.64 

- 6.64 

89.7 

92.3 

-20 

-1.82 

-6.72 

-6.72 

89.7 

92.4 


-15 

-1.82 

-6.70 

-6.70 

89.7 

92.4 


-10 

-1.82 

-6.68 

-6 • 68 

89.7 

92.4 


- 5 

-1.81 

-6.66 

-6.66 

89.7 

92.3 


+ 5 

-1.80 

-6.62 

-6.62 

89.7 

92.3 


* 4 * 1 0 

-1.80 

-6.60 

-6.60 

89.7 

92,2 


+ 15 

-1.79 

-6.58 

-6.58 

89,7 

92.2 


+20 

-1.79 

-6.56 

- 6.56 

89.7 

92.1 


+25 

-1.79 

-6.54 

-6.54 

89.7 

92.1 

X(6) 

0 

-1.81 

-6.64 

-6.64 

89.7 

92.3 


-20 

-1.91 

-6.99 

-6.99 

89.7 

92.7 


-15 

-1.88 

-6.90 

-6.90 

89.7 

92.6 


-10 

-1.86 

-6.81 

-6.81 

89.7 

92.5 


- 5 

-1.83 

-6.72 

-6.72 

89.7 

92.4 


+ 5 

4 1 0 

-1.78 
-1 .76 

-6.56 

-6.48 

-6.56 

-6,48 

89.7 

89.7 

92.2 

92.1 


+ 15 

-1.74 

-6.40 

-6.40 

89.7 

92.0 


* * */ 

+P0 

-1.71 

-6.32 

-6.32 

89.7 

91 .9 


1 4m* 

+25 

-1.69 

-6.25 

-6.25 

89.7 

91.8 


* X ( 2 ) , X(3) and X(4) do not cause appreciable changes in 

the value of u, v, w ,W 1 and to 2 * 
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TABLE 4. OPTIMISATION BESPITS 


Table 4.1. For first starting point 



! 

1 

! Initial value 

1 

.JL* . 

Final value 

Percentage change 

X(1) 

0.5 

0.3656 

21.0 

X(2) 

,•0.1 

0.9743 

2.7 

X(3) 

0.5 

0.4903 

1.4 

3?(X) 

0,5335 

0.4291 

20.0 

rX 1 2 ) 

-0.231 

-0.099 

57.0 

g(l3) 

-0.231 

-0.098 

57.0 

g( H) 

-0.231 

-0.098 

57.0 



Table 4.2. 

For 2nd starting point 



1 

{ Initial value 

i 

t 

1 ! 

{ Final value { Percentage change 

i i 

! _ I 

X(1) 

0.008 

0.005917 

26.0 

X(2) 

0.008 

0.007903 

1.2 

X(3) 

0.050 

0 .048760 

2.4 

F( X) 

0.011503 

.,0.009268 

18.5 

g(12) 

-0.231 

-0.093 

59.0 

g( 13) 

-0.231 

-0.097 

57.0 

g(l4 ) 

-0.231 

-0.095 

58.0 
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APPENDIX A 

The 40 x 40 matrix [_T] ; See Pig. 3 for the definition of 

symbols a, b, c, 9 and 


'1 


l T” 

i ** 


0 

0 

S', 

0 

0 

s, 

0 

0 

s 

0 

•0 


7 


0 

s. 

0 

0 

s, 

0 

0 

s 

0 

0 

s 

0 


7 


0 

0 

s„ 

£ 

0 

0 

S 4 

0 

0 

s < 

0 

0 

s 


8 


M i ! 


i 0 


u 


-b 0 

1 0 
0 1 


b 2 0 0 '0 -xP 0 00 

-2b 0 0 0 3b 2 0 0 0 

0 ~2b -b 0 0 b 2 0 0 


Q 

y> 2ti 


2b 
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1 ~t> 0 b 2 0 o -b 3 0 0 0 

0 1 o -2b 0 0 3b 2 0 0 0 

0 0 10 -b 0 0 b 2 0 0 

-0 0 0 2 0 0 -6b 0 0 0 

0000 10 0 -2b 0 0 

0000 02 0 0 -2b 0 













— 


j 

b 4 

0 

0 

0 

0 


-b 5 

0 

0 

0 

0 

0 

J 

1 

~4b J 

0 

0 

0 

0 


5b* 

0 

0 

0 

0 


j 

1 

0 

-b 3 

0 

0 

0 


0 

0 

0 

0 

0 


[ li 2h\ “ ! 

1 

2 

12b* 

0 

0 

0 

0 

- 

-20b 3 

0 

0 

0 

0 

! 

i 

i 

| 

0 

? 

3b 

0 

0 

0 


0 

0 

0 

0 

0 1 


! 

L 

0 

0 

CM 

CM 

0 

0 


0 

-2b 3 0 

0 

0 


i 

1 

a 

0 

a 2 


0 

0 


a 3 

0 

0 

0 ' 

\ 8-1 = 

0 

1 

0 

2a 


0 

0 


3a 2 

0 

0 

0 


: o 

0 

1 

0 


a 

0 


0 

2 

a 

0 

0 


[ a 4) “ [ S 4a S « ] 
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LV 


[ s s] 


f VI 


■7 ri 


i 

i 

a 

0 

2 

a 

0 

0 

a 3 

0 

0 

— ! 

o i 

1 

0 

1 

0 

2a 

0 

0 

3a 2 

:0 

0 

0 i 

i 

j 

0 

0 

1 

0 

a 

0 

0 

a 2 

0 

5 

0 

! 

f 

t 

0 

0 

0 

2 

0 

0 

6 a 

0 

0 

0 

! 

i 

i 

0 

0 

0 

0 

1 

0 

0 

2a 

0 

0 

! 

\ 

0 

0 

0 

0 

0 

2 

0 

0 

2a 

0 

f 

4 

a 

0 

0 

0 

0 

a 3 


0 

0 

0 0 ; 

f 

X 

4a J 

0 

0 

0 

0 

5a 4 


0 

0 

0 0 













0 

a 

0 

0 

0 

0 


0 

0 

0 0 

12a 

0 

0 

0 

0 

20a 3 


0 

0 

0 0 


0 

p 

3a^ 

0 

0 

0 

0 


0 

0 

0 0 


0 

0 

2a 2 

0 

0 

0 

2a 3 

0 

0 0 

** 

1 

0 

c 

0 

0 

c 2 

0 

0 

0 

° 3 ] 


0 

1 

0 

0 

c 

0 

0 

0 

2 

c 

s 

0 

- 

0 

0 

1 

0 

0 

2c 

0 

0 

0 

_ 2 

3c 

r 

t. 

e* 

u 6a 

s 

6,] 









1 

0 

c 

0 

0 

c 2 

0 

0 

0 

3 1 

C | 

i 


0 

1 

0 

0 

c 

0 

0 

0 

c 2 0 } 

i 


0 

0 

1 

0 

0 

2c 

0 

0 

0 

_ 2 
3c 


0 

0 

0 

2 

0 

0 

0 

2c 

0 

0 


0 

0 

0 

0 

1 

0 

0 

0 

2c 

0 


( 

i o 

0 

0 

0 

0 

2 

0 

0 

0 

6c 



82 


{ S r , 

L 6b J 


! 0 
! o 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

2c 2 

0 

0 


0 

C' 

0 

0 

3c‘ 

0 


c 

0 

4c' 

0 

0 

12c 2 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

2c' 


0 c 

A 


o 

0 


0 4c' 

0 0 


0 

5c 4 

0 

0 

*2 

2c" 


[s 7 ] 


i d-b c (a-b)‘ 
1 * > t > q 


(a-b)c c£ (a-b)^ 


9 


o (a-b)c 
27 ’ 27 


2 ' £ 
27 


27 


I 


8a 


S 8bl 


l 



0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

] - 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 _ 


[ S 8b] 


(5a 4 o, 3a 2 c 3 - 2a 4 c, -2ao 4 + 3a 3 o 2 , 


3 - to 2 o 3 , Sac"* ) 


4 


(5b 


4 3b 2 c 3 _ 2b 4 Cj 2bc 4 _ 5b 3 c 2 ? 


5 - 4b 2 c^, -5bc 4 ) 
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ELEMENTS OP STIFFNESS MATRIX: 

See Eg.. (24) for the definition of the symbols 

m. , n. , p. , q. , r. and s. . 
i’ 4 a’ i i 


k, 




1D l ra j^( m i +m j~2> n i +n j ) + ^/(Pi+PjJ Ii+lj-2) 

+ Hi - f )[ n^E » \ +n j " 2 ) + P ± P ^ (p-j+P j “ 2 , 

q±+q^) ] +| i 1 (1-")'0^P i + ^> q ± ) ] F(m.+p.-1, 
n^+qj+1) + j^(1-W)n jL p^ 1 F(m jL +p..-1 , n ± +q-1 ) 

- ( ) m i I '^ m i +r j~ 1 » n i +s ^ + m j^'( m j +r i “^ » 

n s + 0 i ) ] - ( f« +v V [ < 1 i I ' (p i +r r V 8 ^ + 

^(Pj+rj, 4 j +3 1 - 1 )]- < 1 -’>J) ft; ^ |_ n i p 


q 


riji+s^-l ) + n^+s ± -1 ) + P i P(P i + r j“ 1 » 

q^s.) + P 3 I ‘(Pj +r i“ 1 » q D +S i)] + |_^k + f*Vft + 
2V.f x ^^+ 2(1-P)^] J < r i +r j- S i +S P 

2 ■ 

+ ft Vi ( v 1)(l r 1,f< v r j 4 ’ s i +s j ) 

. t 2 

4- D i G.(s r 1 ) (s r 1 )F(r i+ r., S.+S.-4) + ^ 

[ 2(1 -V) r i r j s i s 3 r i S 3 ( ' r i" 1 ^ S r 1 ^ + ^ )r j s i^ r j" 1 ^ 

(o ± -i 5 '] 



where F(m,n) = c n+1 |^a m+1 . (--b) m+1 
for i & 3 = 1,2,. ..,40. 


.j 


ml nl 
(m+n+2 ) I 


ROTATION MATRIX [R"j: 

1L 


t R j 


0 

0 

6 

o 

o 

o 


0 

R, 

0 

0 

0 

0 

0 


0 

0 

R^ 

0 

0 

0 

0 


0 

0 

0 

R, 

0 

0 

0 


0 

0 

0 

0 

R 1 

0 

0 


0 

0 

0 

0 

0 

R, 

c 

0 


0 

0 

0 

0 

0 

0 

R, 


3 ! 


a. 


cos© 

sin© 


[h] 


0 

0 

ff** 

1 

0 

0 


0 


-c^ 


0 


-c 1 s 1 


0 


C 1 S 1 


0 


- C 1 S 1 


0 
0 
c i 
0 


0 


C 1 S 1 


- S 1 

0 

2 

C 1 

” s 1°1 


0 


C 1 S 1 


0 


C 1 S 1 


0 
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0 0 


0 

0 



0 

0 

0 


2s l c 1 



2 

1 
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appendix b 

MATRIX [a] s See Rig. 4 for definition of various symbols. 
Defining 


b i 

= y 3 - 


c -i 

= X 

X. 


m 

3 

jp 

a i 

= X .V 

- X v . 

1 nr a 

A 

= area 

of triangle. 



a '11 

a 12 

a 13 

a 14 

a 1 5 

a 1 6 

a 21 

a 22 

a 23 

a 24 

a 25 

a 26 


J 


a. 4 »| ***** 

< d i 

+ ~b ± \ + c 



a 1 3 = 

(d 3 

+ b . 1 + c 
D 

!^Vj)/2A 


J3 

U1 

11 

( A, 

+ V- + c 

: m V))/2A 


a 1 2 * 

a 1 1 

a 14 

= a 13 

a 1 6 ~ a 1 5 

a^-j « 

a 22 

»** a.j ^ 

a 23 

= a 24 = a 13 

a 25 “ 

a 26 

= a 15 
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MATRIX 


M= 


1 

1 

~ 2A 

r b . 

1 

0 

b . 

3 

0 

m 

0 

1 0 

c . 

1 

0 

c . 

0 

0 

c 

m 

I 

c . 

1 

b . 

1 

c . 

3 

b . 

3 

c 

m 

b 

m 


MATRIX T03: 


L°] 


1 1 

V 

0 

E 



= 9 ! v 

( 1 -v 2 ) 

1 

0 

( 1 \V\ 

0 

0 

llr-Rl 



2 

EON MATRIX [r]: 



r 

! R i 

0 

0 "j 

[r! = jo 

t 

R 1 

0 

I 0 

L. 

0 

R 1 _ 

j r ii 

1*12 

r 13~ 

f R l"j = , r 21 

r 22 

r 23 

1 r 31 

r 32 

r 33 


= y i 

“ 7 3 5 


z . - z. 

1 0 
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1 = (x? . 

J 31 

+ 4 + ^ 


= . (x 2 . 

v 

2 p 

+ % + *2j) 


'll = X jj/ L 

r 12 = ^i/ 1 

r 13 _ z^/L 

21 = 

r 31 = ^ 

r 32 = Z m/ P 

31 = r 11 r 23 ' 

- r 13 r 22 


32 = r l3 r 21 ' 

mm "Y^ “K' 

11 23 



r 33 = r tl r 22 “ r l2 r 21 



